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ABSTRACT 

Unconfined relativistic outflows from rotating, magnetized compact objects are often well-modeled by as- 
suming the field geometry is approximately a split-monopole at large radii. Earlier work has indicated that 
such an unconfined flow has an inefficient conversion of magnetic energy to kinetic energy. This has led to the 
conclusion that ideal magnetohydrodynamical (MHD) processes fail to explain observations of, e.g., the Crab 
pulsar wind at large radii where energy conversion appears efficient. In addition, as a model for astrophysical 
jets, the monopole field geometry has been abandoned in favor of externally confined jets since the latter ap- 
peared to be generically more efficient jet accelerators. We perform time-dependent axisymmetric relativistic 
MHD simulations in order to find steady state solutions for a wind from a compact object endowed with a 
monopole field geometry. Our simulations follow the outflow for 10 orders of magnitude in distance from the 
compact object, which is large enough to study both the initial "acceleration zone" of the magnetized wind as 
well as the asymptotic "coasting zone." We obtain the surprising result that acceleration is actually efficient in 
the polar region, which develops a jet despite not being confined by an external medium. Our models contain 
jets that have sufficient energy to account for moderately energetic long and short gamma-ray burst (GRB) 
events (~ 10 51 -10 52 erg), collimate into narrow opening angles (opening half-angle 9j w 0.03 rad), become 
matter-dominated at large radii (electromagnetic energy flux per unit matter energy flux a < 1), and move at 
ultrarelativistic Lorentz factors (<yj ~ 200 for our fiducial model). The simulated jets have jjOj ~ 5-15, so they 
are in principle capable of generating "achromatic jet breaks" in GRB afterglow light curves. By defining a 
"causality surface" beyond which the jet cannot communicate with a generalized "magnetic nozzle" near the 
axis of rotation, we obtain approximate analytical solutions for the Lorentz factor that fit the numerical solutions 
well. This allows us to extend our results to monopole wind models with arbitrary magnetization. Overall, our 
results demonstrate that the production of ultrarelativistic jets is a more robust process than previously thought. 
Subject headings: relativity — MHD — gamma rays: bursts — X-rays: bursts — galaxies: jets — accretion, 
accretion disks — black hole physics — methods: numerical, analytical 



1. INTRODUCTION 

Gamma-ray bursts (GRBs), active galactic nuclei (AGN), 
x-ray binaries, and pulsar wind nebulae (PWNe) are among 
the most powerful systems in the Universe. Their power 
originates from a central engine that contains a rotating, 
magne tized compact object suc h as a neutron star or black 
hole dGoldreich & Julian! [19691: IBlandford & Znaiekl [l977h 
or from a surrounding accretion disk ( Shakura & S unyaevl 
119731: iNovikov & Thornd 11973b iLovelacd Il976lh These 
systems obtain their angular momentum and strong mag- 
netic field from their environment either by advection dur- 
ing their formation or through accretion which is known 
to amplify any weak fie ld by magnetorotational turbulence 
(Balbus&Hawley 1998). The region around the compact 
object is often expected to contain a highly-magnetized 
dipolar magnetosphere th at either threads the neutron star 
dGoldreich & Julian! 1 19691) or develops via a c cretion around 
the b l ack hole (IB landford & Znajekl 1 1977b iNarayan etaD 
l2003t iMcKinneyl 120051) . For axisymmetric rapidly ro- 
tating systems, the dipolar magnetosphere can be well- 
modeled by an approximate split-monopole field geometry at 
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large radii once the magnetohydrodynamically-driven (MHD- 
driven) outflow has passed the so-called light cylinder (i.e. 
Alfve n surface) and reaches a point where the flow is uncon- 
fined (Contopoulos et al. 1999; McKinnev 2006b c). For ex- 
ample, astrophysical jets are typically confined by some exter- 
nal medium such as a disk, disk wind, or envelope of matter. 
If the jet remains highly-magnetized far from such confining 
media and passes far beyond the light cylinder, then the mag- 
netic Jie\djgometr£wi\\^^ monopolar 
(e.g. McKinn eyll2006bl:lTchekhovskoy et al.ll2008l) . 

The Crab Pulsar is the quintessential astrophysical ob- 
ject for which the unconfined split-monopole field geome- 
try remains a key model element. One of the most con- 
tentious issues is how to reconcile Crab PWN observa- 
tions with MHD and pair-creation theories. Calculations of 
pair formation fronts both n ear the surface of the neutron 
star i n polar gaps (ISturrocklll971t iRuderman & Sutherland 
19751: f Scharlem ann et alj|1978t iDaughertv & Harding! 1 19821: 
Hibschman & A rons 2001a, b) and farther from the neutron 
star in slot gaps dAronslll983l) and outer g a ps dChengetal 
19761 1198ft iRomani & Y adigaroglu 1995; Ro man! 1 199" 



Cheng et al.l 120001) suggest that the ratio of electromagnetic 



energy flux to matter energy flux in the inner pulsar wind is 
erq ~ 10 3 -10 5 and the Lorentz factor is 70 ~ 10-1000. This 
wind is believed to terminate in a standing reverse shock at 
a distance of ^0.1 pc, i.e., at ~ 3 x 10 11 neutron star radii. 
Observations of the shocked gas, coupled with modeling, in- 
dicate that the pre-shock plasma has a weak magnetization, 
Coo < 0.01, which is 5-7 orders of magnitude smaller than the 



2 



A. Tchekhovskoy, J. C. McKinney, & R. Narayan 



initial m agnetization dRees & Gunn 1974; K ennel & Coronitil 
ll984allbtlEmmering & Chevalieril 1987E 

How does the high-cr wind flowing out of the star con- 
vert essentially all of its Poynting energy flux into ki- 
netic energy flux? This remains an enigma, despite three 
and a half decades of study, and has been coined the 
"cr problem." For the case of a neutron star endowed 
with a split-monopole poloidal magnetic field — a par- 
ticularly simple geometry — it can be shown analytically 
that ideal MHD processes can transfer at most 0.1% of 
the Poynting energy flux from the Crab Pulsar to the mat- 
ter ([Beskin, Kuznetsova, & Rafikov Hl998l) . That is, the wind 
should remain highly magnetized out to the distance of the ter- 
mination shock. In this model, the magnetization near the ter- 
mination shock is expected to be o-qo ~ (ero7o) 2//3 ~ 10 4 3> 1, 
and the Lorent z factor is expected to be 7„o ~ (<7n7n)^ 3 ~ 100 
dMichell fl969l i Goldreich & Julian! [197(3: ICamenzindl [1981 
iBeskin et alJl!998h . This estimate of the magnetization dis- 
agrees with the observationally inferred value of er^ < 0.01, 
and the estimate of the Lorentz factor is far smaller than 
the 7no ~ 10 6 inferred from obser vations dKennel & Coronitil 
1984b; Spi tkovsky & Aronsll200l . 

One might suspect that the above results are artificial, since 
they are derived for the special case of a split-monopole ge- 
ometry. However, an approximately split-monopole is ac- 
tually quite an accurate description of the far regions of a 
dipolar pulsar magnetosphere, and various studies have indi- 
cated that the low efficiency of the split- monopole magneto- 
spher e carries over to the dipolar problem ( Contonoulo s et al.l 
19991; lUzdenskv] [2001 iGruzmovl 12001 iMcKinnevI | 2006c; 
Spitkovsky 20061 IKomissa rov 2006t iBucciantini et al.l 2006). 
This is the reason for continued interest in the split-monopole 
problem. For compactness, we hereafter refer to the case 
where a star is endowed with the split-monopole magnetic 
field geometry as simply the monopole magnetic field geom- 
etry case even though the global solution away from the star 
is not exactly monopolar. 

Various studies have explored the conditions needed 
for strong acceleration of a relativistic magnetized wind 
and effic ient conver si on of magnetic e n ergy to kinetic 
energy . ICamenzindl (11987b. ICamenzindl d 1989b. lLietal.1 
<fl992h . iBegelman & Lil d 19941) . and IChiueh et al.1 d!998l) 
showed that, for efficient energy conversion to occur, 
magnetic field lines should expand away from one another 
and away from the equatorial plane. This field geom- 
etry was identified as a "magnetic nozzle" because the 
expansion of field lines away from the equatorial plane 
is geometrically similar to the expanding outer edge of 
nozzles (e.g. de Laval nozzle) intended to launch a su- 
personic flow. By this argument, the monopole geometry 
is part icularly inefficient since field lines are perfectly ra- 
dial dBeskin et al.ll 1998b lBogovalovll2001l lKomissarovll2006b 
Bucciantini. Thompson. Arons. Ouataert. & Del Zannal 
2006b Barkov & KomissarovT 2008; 

Komissarov, Vlahakis, Konigl, & Barkov 2009). Field 
geometries other than monopolar/dipolar do manage to 
convert Poynting flux to kinetic energy flux more effi- 
ciently, reaching a < 1 dLT. Chiueh. & Begelmanl 119921: 
I Begelman & Lil 1 1994b I Vlahakis & Koniglll2003allbh IVl ahakis 
2004b lBeskin& Nokhrina l2006i iBarkov & Komissarovl2008b 
IKomissarov et al.l 120091) . However, all the cases considered 
so far that show efficient acceleration to large Lorentz factors 
(7 > 100), have involved outflows that were restricted to flow 



inside coUimating walls with prescribed shapes or confining 
pressure profiles that induce collimation. Some prior ideal 
MHD work claiming to solve the er-problem prescribed 
the field line shape, i.e., did not have a self-consistent 
(global force-balanced) solutio n dTakahashi & ShibataiEHH 
IContopoulos & Kazanasl 120021: iFendt & Ouyedl 12004 . It 
remains unknown whether these "jet" models continue to 
exhibit efficient acceleration if the walls are either removed 
or given a different shape or if the field line shape is 
self-consistently computed. 

By consi dering small pertur bations to the monopole field 
geometry, Beskin et al. (1998) derived self-consistent solu- 
tions of highly magnetized monopole o utflows near the mid- 
plane a nd found inefficient acceleration. Lvubarsk v & Eichlerl 
d2001l) have extended their analysis to the polar regions of 
such outflows and showed that highly-magnetized (a 3> 1) 
collimated relativistic jets can form there, however, they did 
not explore whether these jet s can become mat te r-dominated 
a < 1 ) at relevant distanc e s. lBogoyalovld2001l); IK omissarov 
2006); IBucciantini et al.l d2006l) : IKomissarov et al.l (12009b 
have numerically simulated unconfined magnetized out- 
flows and confirmed their low efficiency at converting 
magnetic to kinetic energy, i.e. the outflows remained 
highly magnetized g > 1 ou t to the simulated distances. 
Tomimatsu & Takahashi (2003) found solutions to cold ideal 
MHD jets that were limited to lie inside very narrow 
boundaries with 9j <C 1/co a nd in which the poloidal cur- 
vature force was neglected dVlahakisI 120041) . Recently 
Zakamska, Beg elman. & B landfordl (I2008h studied conver- 
sion of internal energy to kinetic energy in hot ideal MHD 
jets, assuming a purely toroidal magnetic field and also 
assuming self-similarity that does not allow for efficient 
conversion of magnetic to kinetic energy. We note that 
studies of highly magnetized flows in the force-free ap- 
proximation (whi ch n eglects matter inertia and kinetic en- 
ergy, § §13.11 and 15. Il l have given much insig ht into how 
jets/winds are launched and into their structure dCamenzind] 
| l987t lAppl & Camenzindlll993b IContopouloslll995b|: i Fendtl 
1997t lLovelace & Romano val |2003b lLovelace et all l2006t 
Uzdenskv & MacFadyenll2006L 120071: iMcKinney & Narayanl 
l2007bt iNarayan et aljl2007t iTchekhovskoy et al.ll2008l) . 

In addition to the above studies, various models have been 
propo sed that involve dissipative processes, e.g. reconnec- 
tion dLvutikov & Uzdenskvll2003b IZenitani & Hosh1nol l2008b 
Malvshkin 2008; Uzdenskv 2009), as possi ble resolutions 
to the cr-problem. In a stri ped wind model ( Micheil fl 97 lb 
ICoronitil[l9 90; Michel Il994b reconnection in a warped equa- 
torial current sheet converts magnetic energy into the ki- 
netic energy of the plasma. However, it remains uncer- 
tain whether such a reconnection process is fast enough to 
accelerate the plasma as required dLvu barskv & KirkJ 12001b 
Kirk & Skiseraasen 2003; iLyubarskyl 12005b . More recently, 
Petri & Lvubars kv (2007) have shown that magnetic recon- 
nection of the warped equatorial current sheet may occur 
right at the wind termination shock, leading to a decrease in 
the inferred pre-shock wind magnetization. Bege lmanl d 1998b 
suggested that toroidal field instabilities may lead to dissipa- 
tio n of the toroidal fiel d, th us circumventing the argu ments 
by iRees & Guniil d 1974b and iKennel & Coronitil dl984al) that 
inferred a low value of a in the termination shock by as- 
suming the shock contains an ordered and purely toroidal 
field. New analysis is required to check consistency be- 
tween observations and MHD models involving disordered 
toroidal fields due to MHD instabilities. However, even if 
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linear MHD instabilities are present, as argued by Begelman 
(1998), their effectiveness remains unknown since they may 
evolve to a saturated non-linear state that has negligible dissi- 
pation, e. g. as demonstrated recen t ly for outflows from black 
holes bvMcKinnev & Blandfordl (120091) . We note also that 
iNarayan. Li. & Tchekhovskoy (2009) showed for a simple jet 
configuration that the linear instability growth rate is much 
lower than one might expect from standard instability criteria. 

In this paper, we present a detailed study of the relativistic 
magnetized monopole wind using both numerical and analyt- 
ical ideal relativistic MHD methods. The simulations we re- 
port here involve a much larger dynamic range than any pre- 
vious published work; they extend in radius from r = 1, the 
surface of the neutron star, to r = 10 10 . The wide range of 
radius allows us to study the solution far into the asymptotic 
region of the wind where acceleration has practically ceased. 
Also, we consider a number of different prescriptions for the 
mass-loading of field lines at the stellar surface. 

The goal of this study is two-fold. First, we wish to focus on 
field lines in the equatorial region of the outflow to study the 
classic a problem. In particular, we wish to compare numeri- 
cal results with previously published analytical results for the 
asymptotic Lorent z factor and magnetization parameter (e.g., 
Beski net al.lll998l) . Second, we wish to study the behavior of 
field lines near the rotation axis. Even though the monopole 
problem is highly idealized, nevertheless, we believe the polar 
field lines in this model may be viewed as analogs of relativis- 
tic jets and indeed may even be directly relevant to relativis- 
tic jets that become unconfined at large radii. Our goal is to 
understand if there are any limitations on acceleration along 
polar field lines. In other words, is there a a problem for jets? 

In $2] we describe the problem setup and the numerical 
method we use to carry out the simulations. In $3] we present 
results for two simulated models: M90 and M10. In $4] we 
study the shapes of field lines and explore the connection be- 
tween field line shape and acceleration. We show that there is 
a large difference between equatorial and polar field lines. In 
$5] we investigate what role if any is played by signals travel- 
ing from one region of the magnetosphere to another, and how 
this affects the efficiency of acceleration. Once again, we dis- 
cover that equatorial and polar field lines have qualitatively 
different rates of acceleration and efficiency. We discuss the 
implications of our results in §|6] and conclude in ^7] 

We work throughout with Heaviside-Lorentz units, and we 
set the speed of light, the radius of the central compact ob- 
ject and the radial component of the surface magnetic field to 
unity. We use spherical polar coordinates, r, 9, tp, as well as 
cylindrical coordinates, R = rsind, ip,z = rcosO. 

2. PROBLEM SETUP 

2.1. Initial Conditions and Time Evolution 

We idealize the central neutron star as a perfectly conduct- 
ing sphere that we refer to as the "star." We assume that the 
star has a split-monopole magnetic field configuration, with 
unit field strength at the stellar surface (this choice sets the 
energy scale). Exterior to the star, we initialize the system 
with a low but finite rest-mass density atmosphere, which is 
done because the code cannot accurately evolve a large den- 
sity contrast between the initial and injected density or a large 
a value in the initial atmosphere. The density of the atmo- 
sphere is chosen so that it is dynamically unimportant (the 
kinetic energy of the piled-up atmosphere is much less than 
the kinetic energy of the wind). This atmosphere is easily 



swept away by the outflowing MHD wind and has no effect 
on the final results. This was confirmed by considering other- 
wise identical models but where the atmosphere density was 
20 times lower. We find all our results are converged indicat- 
ing negligible impact by the atmosphere on the injected wind. 

The initial system has no rotation, so field lines are perfectly 
radial and both Bg and B v vanish. Starting with this initial 
configuration, we impose a uniform rotation on the star and 
study the time evolution of the external magnetosphere. As 
the star spins up, the footpoints of magnetic field lines are 
forced to rotate, and this generates a set of outgoing waves 
traveling at nearly the speed of light. A short distance behind 
the outgoing wavefront, the magnetosphere settles down to a 
steady state. We are interested in the properties of this steady 
MHD wind. 

The computational domain in our simulations is the upper 
hemisphere, < 9 < n/2, with radius extending from the 
surface of the star, r = 1, to an outer edge at r = 10 12 . We 
note that most calculations in the literature are limited to a 
very small radial range due to the need to always resolve the 
time-d ependent compact object (e.g. iMcKinnev & Naravanl 
l2007al) . Our choice of a very large range of radius allows us 
to study both the initial "acceleration zone" of the magnetized 
wind as well as the asymptotic "coasting zone." 

2.2. Boundary conditions 

At the polar axis, 9 = 0, and at the midplane, 9 = ir/2, we 
use the usual antisymmetric boundary conditions. At the outer 
radial boundary (r = 10 12 ) we apply an outflow condition. At 
the stellar surface, we set the poloidal component of the 3- 
velocity of the wind to a fixed value v p directed along the 

poloidal magnetic field B p . We choose v p = 1/2 in all the sim- 
ulations reported here. Also, we choose the angular velocity 
of the star to be il = 3/4 (i.e. v v = 3/4 at the stellar equator). 
These boundary conditions are the same for all simulations. 

We assume that the magnetic field is frozen into the star. 
The radial component of the field is continuous across the 
stellar surface, so we enforce the boundary condition B r = 1 at 
r = 1 . Since we inject a sub-Alfvenic flow at the surface of the 
rotating star, Alfven and fast magnetosonic waves communi- 
cate from the magnetosphere back to the surface and generate 
self-consistent non-zero values of Bg and B^. The fluxes at 
the stellar surface are set by using "outflow"-type boundary 
condition (i.e., flowing out of the computational domain into 
the star) on Bg and B v . 

Since we fix v p and fl at the stellar surface, the toroidal 
component of the 3-velocity v v is determined by the condition 
of stationarity: 

v tp = m+B tp v p /B p . (1) 

This equation follows by decomposing the wind velocity into 
rotation with the field line (the first term) plus motion parallel 
to the field line (the second term). 

The final boundary condition at the stellar surface is the 
plasma density. This controls how much mass is loaded onto 
field lines at their footpoints. The different simulations we 
report in this paper correspond to different choices for p(9f p ), 
the profile of density as a function of polar angle 6*f p of field 
line footpoints across the stellar surface. 

2.3. Numerical Approach 

We solve the time-dependent axisymmetric equations of 
special relativistic MHD (ignoring gravity) at zero tem- 
perature with second order accuracy. For this we use the 
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TABLE 1 
Simulation Parameters 



Name 




/^max 


fi ollt r Q f c n S[°] x 2 


Resolution 


Eft'. Resolution 








Constant density-on-the-star models 






M90 


90 


460 


— 5 0.25 4 11.25 1/4 


1536x384 


3918x768 


1000M90 


90 


1000 


— 100 0.25 3 90 1 


2048x384 


2784x384 








Variable density-on-the-star models 






M45 


45 


460 


40 100 0.25 3 90 1 


2048x384 


2784x384 


M20 


20 


460 


40 0.55 100 0.05 3 20 1/3 


3072x384 


3216x576 


MIO 


10 


460 


40 5 0.25 4 10 1/3 


3072x384 


7837x1152 






Variable density-on-the-star + wall at 8 = 6> m ax 




W10 


10 


460 


40 0.7 100 0.25 3 90 1 


2048x128 


2786x128 


W5 


5 


460 


40 0.55 100 0.25 3 90 1 


6144x128 


8358x128 



code HARM dGammie et all 120031: iMcKinnev & Gammid 
120041: livIcKinnevi l2006bl) w ith recent improv ements 
(Mign one&McKinneyl 120071: iTchekhovskoy et all 120071 
120081) . We use HARM's second-order MC limiter method 
for spatial interpolations and a second-order Runge-Kutta 
time-integration. To improve the accuracy of the simulation, 
before each reconstruction step we interpolate the ratio 
of the numerical to approximate analytical solution, as 
described in Appendix [C] To speed up the computations, 
we stop evolving regions of the s olution where the wind 
has achieved a steady s tate (c.f. Komiss arov et al.l 120071 
ITchekhovskoy et aT] 120081) . This technique allows a large 
gain in speed of up to a factor ~ 10 10 , proportional to the 
maximum radius of the simulation. Also, since our interest 
is in cold flows in which the plasma internal energy and 
pressure are negligibly small, we set these quantities (and all 
their derivatives) identically to zero and ignore the energy 
evolution equation. For example, the conversion of conserved 
to pri mitive quantities is performed identically to the thermal 
case dMignone & McKinnev 2007), but pressure and internal 
energy and all their derivatives are set to zero. 

HARM is a flexible code that permits the use of an arbitrary 
coordinate system. We employ a radial grid in which the res- 
olution is very good near the star, but the cells become much 
more widely spaced at large radii. In terms of a uniformly- 
spaced internal code coordinate x\, we write the radial coor- 
dinate of the grid as 

r(x\) = ro + exp \x\ +cH(x\ —x\) x (x\ — x{y\ , (2) 

where H(x) is the Heaviside step function. The lower (up- 
per) x\ value is determined by the lower (upper) radial edge 
of the grid: r = 1 (r = 10 12 ). The parameters ro, r = r(x\), 
c > 0, and n > 2 allow us to control how resolution varies 
with r. For instance, at ro <C r < r the grid is logarithmic, 
with near-uniform relative grid cell size Ar/r k, const. At 
r = r the grid smoothly switches to hyper-logarithmic, with 
Ar/r oc In 1-1 '" r for r ^> r. Using a hyper-logarithmic grid at 
large radii is sufficient since all non-trivial changes in quanti- 
ties (e.g. the Lorentz factor and collimation angle) occur log- 
arithmically slowly with r. 

For the angular coordinate, we map 6 to a uniform code 
coordinate X2, which goes from to 1, according to 

e(x 2 )=x 2 l + H(x 2 -x 2 ) X -f$S ■ 0) 

x 2 \2 x 2 J \l-x 2 J 

In this grid, a fraction x 2 (typically ~ 1 /2) of the cells are dis- 



tributed uniformly between 9 = and 9 = 9 and the remaining 
cells are distributed non-uniformly at larger angles. 

In different models we utilize grids with different values of 
the parameters ro, r, c, n, x 2 , and 9. Table [Ugives the details 
and estimates the effective resolution of our models in terms 
of the typically-used uniform angular grid and exponential ra- 
dial grid (i.e. r(x\) = ro + exp[xi]). 

3. SIMULATION RESULTS 

3.1. Baseline Model M90 

For our fiducial baseline model, called M90, we set the den- 
sity p(9f p ) at the surface of the star to a constant value equal to 
0.001914, independent of 6, such that B 2 r /(p c 2 ) « 522. The 
rotation of the star launches a relativistic magnetized wind 
that accelerates the input mass. The simulation is run for a 
time 10 12 in units of the light-crossing time of the star. By 
this time, the wind reaches a steady-state solution that is in- 
dependent of initial transients out to a radius 10 10 , so all our 
results are reported out to this radius. The top two panels in 
Fig. Q] show the results. 

Consider first the shapes of the field lines in the poloidal 
plane. We see that, over most of the solution, the field 
lines are only slightly perturbed from their initial purely ra- 
dial configuration. This might be surprising since rotation 
causes the toroidal component of the field to grow substan- 
tially. In fact, over most of the solution, B v is tens of times 

larger than the poloidal field B p = B 2 r +B 2 g , and one might 

think that the hoop stress of the strong toroidal field would 
cause substantial collimation of the field lines. This does 
not happen because the relativistic wind sets up an electric 
field, and an associated electric force per unit volume p e E 
(charge density p e ), that almost exactly cancels the mag- 
netic hoop stress. This is a unique feature of relativistic 
win ds. It has been explored in d etail bv lNaravan et alj (12007b 
and lTchekhovskov et al.l {2008, hereafter TMN08) in connec- 
tion with winds in the force-free app roximation dMcKinnevI 
I2006ak IMcKinnev & Naravanll2007bl) . 

While it is true that the distortion of field lines in the 
poloidal plane is small, nevertheless, there is some distortion. 
It is especially obvious near the rotation axis, where we see 
field lines converging towards the pole and bunching up. In 
fact, very close to the rotation axis, the distortion appears to 
be quite large. 4 Thus, the results indicate that the polar re- 

4 Note that Fig.fT]uses a logarithmic radial coordinate and thus exaggerates 
the effect. The actual change in 8 of a field line is not very severe even at very 
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FIG. 1. — Results for models M90 (upper panels) and M10 (lower panels). Radial dark lines show, for the steady state solution, the poloidal structure of field 
lines from log 10 r = (the surface of the star) to log 10 r = 10. The white radial lines on the outside show the initial non-rotating monopolar configurati on o f field 
lines. The three thick solid lines correspond to the Alfven surface (innermost), the fast surface (middle), and the causality surface (outermost, see j|5.2t . The 
fast surface and the causality surface touch each other at the midplane (8 = tt/2). Top-left: color-coded Lorentz factor 7 in model M90. Top-right: color-coded 
energy flux per unit mass flux fi in M90. Bottom-left: 7 in M10. Bottom-right: fi in M10. In the figure we use "logarithmic" spherical polar coordinates 
(r; = 1 + log 10 r ,81 = 8). The numbers along the horizontal (vertical) axis correspond to the values of log| r along that axis, i.e., r; cos0; — 1 (r; sinf?/ — 1). Even 
though in these "logarithmic" polar coordinates the flow appears to over-collimate toward the axis (e.g., the two leftmost field lines on the upper panels), in fact 
dR/dz remains positive everywhere in the solution. We note that the shape of the shown portion of the Alfven surface is very close to a cylinder R = 1 /O = 4/3. 
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gions of a rotating monopole wind are qualitatively different 
from the equatorial regions of the outflow. 

We next consider the acceleration of the wind. The top left 
panel of Fig.Q]shows the Lorentz factor 7 as a function of po- 
sition. We see that there is an inner acceleration zone extend- 
ing out to a radius r ~ 10 2 , where the Lorentz factor increases 
from its initial small value (« 1.2 at the stellar surface). Be- 
yond this we find a large coasting zone where there is very 
little change in 7. Over most of the coasting zone, the flow 
reaches ultrarelativistic Lorentz factors of 7 ~ 40; however, 
this value is much less than expected if all the available free 
magnetic energy were to be used. 

An axisymmetric magnetized wind has several conserved 
quantities along field lines. Two of these are the enclosed 
magnetic flux $ and the angular velocity fi. A nother is the ra- 
tio of poloidal magnetic flux to rest-mass flux (|C handrasekhar 
[T95llMestell[T96ll: iLi et al.lfl99l lBeskinlll997l) : 

7y($) = 1£0l = const, along field line. (4) 
B p 

Yet another conserved quantity is the quantity p, which 
is the ratio of the tot al energy flux to the rest-mass flux 
dChandrasekharj [1956. iMestell Tl96lt [Lovelace et all [19861: 



IBegelman & Lilll994t[Beskinlll997HChiueh et al.1 11998): 



S + IC E\B v \ + j 2 pv p 



K 



IPVp 



■ const, along field line. (5) 



Here, <S is the Poynting flux, JC is the mass energy flux (rest 
mass 1Z plus kinetic energy), B v is the toroidal field, 



E = Q,RB n 



(6) 



is the poloidal electric field, and p is the mass density in the 
comoving frame of the fluid. The denominator of equation 
(O is the rest mass flux. Since we consider highly magnetized 
winds, we have S ^> K, at th e surface of the star. Moreover, at 
r = 1 in the monopolar flow dMichellll 9691 1 19731) . 



\B V \ ss£ « fisinflfp, 



(7) 



where 6f p represents the value of 9 at the footpoint of a field 
line. Also, at the footpoint, v p and p are constant, and 7 « 1. 
We thus have 5 



p(9{ p ) ex sin 



7fp- 



(8) 



for // 3> 70 (70 is the initial Lorentz factor at r = 1), i.e., p is a 
rapidly increasing function of 6f p . This can be seen in the top 
right panel in Fig. Q] and also in Fig. [2] 

As mentioned above, the quantity p is conserved along each 
field line. However, the two energy contributions to p, the 
Poynting flux S and the mass energy flux K,, are not individu- 
ally conserved. In fact, the outflowing wind converts Poynting 
flux to mass energy, thereby accelerating the wind and caus- 
ing 7 to increase. 

If the wind were maximally efficient at accelerating the 
matter, we would expect S — > at large distance from the star. 
The Lorentz factor would then achieve its maximum value 



7n 



■p. 



(9) 



In practice, the wind falls far short of this maximum. As 
Figs. [1] and [2] show, p is quite large in model M90, with a 
value of 460 at the equator (6>f p = tt/2). However, the Lorentz 

5 In this paper we label field lines by either the enclosed magnetic flux $ 
or the polar angle at the footpoint 8fp. 



factor of the wind, even at a radius of 10 10 , does not exceed 
40. Thus, a magnetized monopole wind is very inefficient at 
accelerating the gas. For an equatorial field line in M90, the 
efficiency factor 7//! is only about 0.09, as shown in Fig. [2] 

Another way of describing the efficiency of conversion 
of energy from electromagnetic to kinetic form is via the 
magnetization parameter a, which is the ratio of Poynt- 
ing to mass energy flux {kennel & Cor oniti 198 4at ILi et al 
ll992H Begelman & Lill994HVlahakis 2004; Kom issarov et al 
1200% . 

a-i = m. (10) 



K 



l 2 pv P ' 



Substituting into equation (|5]l we see that the conserved quan- 
tity p is related to a by 



p = 1 (a+X). 



(11) 



The smallest value possible for a is zero. Therefore, the max- 
imum value of 7 is 7 max = p (eq. |9). 

The magnetization a is not conserved along a field line. At 
the surface of the star, where 7 w 1 , we have <r ~ p — 1 . As 
the magnetized wind flows out and energy is transferred from 
Poynting to matter energy, 7 increases and a decreases. An 
efficient wind would be one in which a asymptotes to a value 
< 1, so that the outflowing material is able to convert at least 
half of its energy flux into matter energy. The numerical so- 
lutions shown in Figs. Q]and[2]fail t0 satisfy this criterion by 
a large factor in the equatorial regions. This implies there 
is no ideal, axisymmetric MHD solution to t he er-problem 
dRees & Gunrill974lKennel & Coronitill 1 9 84allbh . 

There is, however, one promising feature in the results: the 
polar regions of the wind are efficient, with j/p — ► 1 and 
a <C 1 at large distance from the star (Fig. [2j. This interesting 
feature of the monopole problem has not been emphasized in 
the literature. Most previous analyses and discussions have 
focused on equatorial field lines where the efficiency is, in- 
deed, too low to solve the a problem. 

Unfortunately, the actual Lorentz factor along polar field 
lines in M90 is only ~ 20 since this is the value of p for these 
lines. Would we continue to have high efficiency in the polar 
region even with larger values of pi In particular, is it possible 
to have acceleration with high efficiency up to Lorentz factors 
7 > 100, as observed for instance in gamma-ray bursts? For 
this we need to study a model with larger values of p near the 
pole. We describe such a model in the next subsection. 

3.2. Model M10 

From equation (0 we see that an obvious way to increase 
p is to lower the density of the wind at the stellar surface. For 
instance, if we were to reduce p by a factor of ~ 30 relative to 
M90, then we would have a model with p ~ few hundred for 
a field line with 6f p ~ 10° (see Fig. [3). We could then explore 
acceleration along this field line and determine whether or not 
the outflowing wind achieves a coasting 7 > 100. 

Since p varies as sin 2 6f p , this approach would lead to ex- 
tremely large values of p at the equator. As a result, the model 
would require very large resolution to simulate accurately and 
would be extremely expensive. Therefore, for numerical con- 
venience, we consider a model in which we choose the profile 
p(f?f p ) such that p is large near the pole, reaches a maximum 
/"max at a specified footpoint angle f?f p = (9 max , and then de- 
creases with increasing 9f p to an outer value p out at 6f p = ir/2. 6 

6 Decreasing near the equator is actually a physically reasonable ap- 
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FIG. 2. — The panels on the left correspond to model M90 and those on the right to M10. The solid black lines in the top panels show the scaled total energy 
flux ji as a function of sinf?f p , where #f p is the polar angle corresponding to footpoints of field lines. The five colored lines show profiles of 7 at five different 
distances from the star: r = 10 2 (magenta, dash-dotted), 10 3 (blue, long-short-dashed), 10 4 (green, short-dashed), 10 6 (orange, long-dashed), 10 9 (red, solid). 
The closer 7 gets to fj, the more efficient the acceleration. The middle panels show the acceleration efficiency 7//1 and the lower panels show the magnetization 
parameter a = (/i/7)— 1 at the same five d istances. Note that ■y/ fi is largest and a is smallest closest to the pole (sin£>f p C 1). The various dotted lines correspond 
to the analytical model described in j)4.1l The shaded area in the bottom panel indicates the regions of the solution where mass energy flux exceeds Poynting flux 
(it < 1, or f/ fi > 1 /2). Note that magnetic flux gradually converges toward the polar axis as a function of increasing radius, so the lines corresponding to larger 
radii truncate at small values of 9{ p . 



To achieve this, we choose the density profile on the star to be 

p(#fp) = A) + Pisin Q f? fp , (12) 

where the constants pg, p\ and a are adjusted so that the model 
has the desired values of p m - dx , 9 m . dx and p out . 

We have simulated a series of such models, all with /i max = 
460 and p out = 40 7 , and with different values of # max , viz., 
45°, 20°, 10°, which we refer to as models M45, M20, M10, 
respectively. Figure [3] shows the profiles of p(9f p ) and p(6f p ) 
for the model M10. Parameters of the various models are 

proach to model the equatorial region if it were to contain a weakly magne- 
tized pulsar current sheet or an accretion disk. 

7 We have confirmed that the precise value we choose for fj, ou t 2> 1 is 
unimportant so long as we are only interested in acceleration along polar 
field lines. 



summarized in Table Q] The models are well-converged, with 
field-line invariants f2($), r/(3>), and p($) conserved along 
field lines to better than 15%, as Fig.|4]shows for, e.g., mod- 
els M90 and M10. 

The bottom two panels in Fig.[T]show results corresponding 
to M10. As before, we see that the poloidal structure of the 
field is largely unaffected by rotation. There is of course some 
lateral shift of field lines, the effect being larger near the pole 
than near the equator, with maximum field line bunching close 
to the axis. 

The Lorentz factor distribution again confirms the trends 
seen in M90. Near the equator, the asymptotic 7 is only ^12, 
giving an inefficient flow with j/p ~ 0.3. The efficiency in- 
creases near the axis and becomes practically equal to unity 
close to the pole; equivalently, a becomes much less than 
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FIG. 3. — The red lines in the two panels show the distribution of p and p as 
a function of footpoint angle 9f p in the baseline model M90. The blue dotted 
lines show the corresponding profiles for a hypothetical model in which p is 
decreased by a constant factor, leading to an increase in p by the same factor. 
Since this causes p to become very large at the equator, leading to numerical 
difficulties, we consider the model M10 (dashed blue lines) that has the same 
behavior of p and p near the axis but is less extreme near the equator. In 
both models M10 and M90 the matter-dominated part of the flow p ~ 1 at 
the surface of the star is resolved on the grid by at least several grid cells. 

unity for these field lines. Most interestingly, Lorentz fac- 
tors nearly as large as ~ 200 are obtained for these field lines. 
In other words, there is no a problem near the axis and it is 
possible to obtain quite large Lorentz factors in this region of 
the outflow. 

3.3. Summary of Key Results 

From the results shown in Figs. [T] [2] we conclude the fol- 
lowing: 

1. Field lines near the equator of a rotating monopole 
largely retain their monopolar configuration out to large 
radii, whereas lines near the pole tend to bunch up 
around the axis. 

2. The acceleration efficiency j/fj, of a rotating monopole 
magnetosphere is low (<c 1) for field lines in the equa- 
torial region, but quite high (~ 1) for field lines near the 
pole. 

We would like to develop an understanding of the physics be- 
hind of these effects. We would also like to know how the two 
effects are related to each other. This is the topic of the next 
two sections. 

4. FIELD LINE BUNCHING 

4.1. Relation to Acceleration Efficiency 

As plasma streams along field lines, relativistic effects be- 
come important near the so-called light cylinder, R = R^ = 
l/f2, where the co-rotation velocity CtR equals the speed of 
light and E = B p (see eq. [6]). As we show in Appendix lAl 



far outside the light cylinder, where VlR ^> 1 and 7 ^> 70, the 
plasm a simply drifts perpendicular to E and B a t the drift ve- 
locity (iBeskin et al.ll 1998.12001 IViahakis 20Ql, 



ExB 



B 2 



E 
B' 



(13) 



The corresponding Lorentz factor is 

2 2 _ B 2 
7 ~ 7dr ~ B 2 -E 2 ' 

Near the star this formula becomes inaccurate since the 
plasma moves at the initial Lorentz factor, 



(14) 



7 



(15) 



In Appendix |A] we show that a combination of these two for- 
mulae does a very good job of describing the Lorentz factor 
at all distances from the star: 



7 2 = 7o-l+7dr 



(16) 

In the asymptotic region of a relativistic outflow (where 
flR ;§> 1,7> 1), we have according to ( fl4] i. 



\B U 



(17) 



Further, this relation is also true at the surface of the central 
compact star (for the monopolar flow, see eq. 0. This al- 
lows us to write the difference between the maximum allowed 
Lorentz factor 7 max = /i and the local Lorentz factor 7 in the 
following convenient form (the numbers in parentheses refer 
to the equations used to derive this result): 



M"7 



E,CD,Q3 E 2 @} (flRB p ) 2 



IPVp 



" 2($) « R 2 



(18) 



By dividing this equation by itself as evaluated at the foot- 
noint and annroximatinp -v— T(_ sa 7 and fJ, — 7f p ~ fi, 



point, and approximating 7 
tain 



2 
fi 



-7f P 
B P R 2 



[B P R 2 ]i v 



1-A 

flfp 



, we ob- 



(19) 



where we have used the subscript "fp" to denote quantities 
evaluated at the field line footpoint 8 and have defined the 
quantity 

a = B p R 2 . (20) 

Equation (fT9l demonstrates that, in order to convert an ap- 
preciable fraction of the total energy flux /1 along a field line 
into matter energy flux (7 times the mass flux), the quan- 
tity a has to decrease appreciably from its initial value at the 
footpoint. This result is known ( c.f. iBegelman & Lil Il994t 
IChiueh et al]|1998tlVlahakisll2004l) . though we have not seen 
as simple a derivation as the one given above. 

For a precisely monopole field, B p oc l/R 2 along each field 
line. Therefore, a = B p R 2 is constant along a field line and 
so no efficient acceleration is possible. This explains why 
acceleration is so difficult in the monopole problem. In order 
to permit acceleration, field lines must move in a cooperative 
fashion transverse to one another so as to allow a to reduce 

8 We note that for field geometries other than monopole, equation (T7) in 
general b reaks down at field line footpoints but holds at the fast magnetosonic 
surface ( j|5.U . Due to this reason, for field geometries other than monopolar 
the subscript "fp" indicates quantities as measured at the fast magnetosonic 
point. 
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FIG. 4. — The panels on the left correspond to model M90 and those on the right to M10. The five colored lines correspond to r = 10", 10 , 10 , 10 , and 10 
(see Fig.[3]for a description of line colors and types) and show profiles of quantities that are preserved along the field lines: n(9f v ) (the top panel), rj(8f p ) (the 
middle panel), and f!(0f p ) (the bottom panel). These lines nearly overlap each other, as expected for a flow in steady state. In the numerical simulations, each of 
the quantities is preserved along field lines to better than 15%. A few grid cells from the polar axis, within which unresol ved m agnetic flux accumulates, are least 
accurate (not shown); but as we have checked, this does not affect the quality of the solution at larger angles (see end of 34.41 . 



with increasing distance from the star. We discuss how this is 
accomplished in the next subsection. 

Meanwhile, as an aside, we describe here an improvement 
to the approximate result $1% which gives the correct numeri- 
cal value of efficiency in a precisely monopolar field. For such 
a field, the Lorentz factor at as ymptotically large distances is 
(iMichelll 1 9691: ICamenzindll 1 9861) 



7 



radial = ^1/3 <<; 



This gives an extremely inefficient acceleration, 



-2/3 



(21) 

7/M ~ 



<C 1. However, the efficiency is not zero as equa- 



tion ( fl9b might suggest, so we need a more accurate version 
of dTTb . According to (fl4"l >. 



QRB„ 1 + 



2j : 



(22) 



Using this equation, assuming B p R = const, along a field 



line, and introducing 7 

R^oo 

~ loo + 



■ lim R - 



we obtain in the limit 



7j($) 



BpR 1 



1 + - 



27? 



(23) 



Now, assuming that the system settles down to a state with 
minimum total energy flux /i($) (equivalent to the mini- 



mal torque condition of Michel [19691) . we find the terminal 
Lorentz factor that minimizes the right hand side of this 
equation: 



, radial 



To 



ft 2 ($) 
r](<f>) 



1 



BpR 1 



1/3 



(24) 



where the approximate equality comes from equation < TT~8T > 
evaluated at the footpoint. This reproduces (f2TT> and shows 
that indeed, for a precisely radial flow, only a small (but still 
non-zero) fraction of the Poynting flux is converted to the ki- 
netic energy of the matter. This derivation was performed for 
precisely monopolar field lines. In actuality the shape of the 
poloidal field lines is slightly changed from radial. The effect 
is small for equatorial field lines and the order of magnitude 
estimate ( 1241 continues to hold. The deviations are larger for 
polar field lines and the Lorentz factor obtained along these 
lines is very different from (F24l >. 

4.2. Field Lines Near the Midplane 

Let us first apply formula ( fT9l to a field line near the mid- 
plane. For acceleration to be efficient, a must decrease along 
the field line, i.e., B p must decrease faster than 1 /R 2 . This can 
be accomplished by moving field lines away from the equa- 
tor towards the axis. Consider a field line with its footpoint 
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located at a small angle 6L = 7r/2-0f p <C 1 from the equator. 
For this field line we can write (T% as 



7 K 



/' 



l-Jt 

m > 



(25) 



where 6' is the polar angle of the field line at a large distance 
from the star, where the Lorentz factor is 7. Here we used the 
fact that $' w 2na9' w 2iraf p 0' fp , where $' is the amount of 
flux enclosed between the field line and the midplane. 

For a nearly monopolar configuration of the field in which 
9' k, 9' fp , obviously we will not have much acceleration 
(7//1 <C 1)- In order to obtain high acceleration efficiency in 
the equatorial region, field lines must diverge from the equa- 
tor so that the 9' values of field lines increase with distance. 
For this to happen, the rest of the magnetosphere must col- 
lectively move away from the equator towards the pole. For 
reasons that are discussed in §0 this does not happen, and so 
acceleration efficiency is at best modest near the equator. 

4.3. Polar Field Lines 

The story is quite different for field lines close to the axis 
(9 <C 1). Consider the initial undistorted monopole configura- 
tion. At the surface of the star, B p is constant, and [B p R 2 ]f p is 
simply equal to $ /ir, where $ is the flux interior to the field 
line. For a pure monopole, this relation is valid at any dis- 
tance from the star, i.e., B p R 2 = <£>/7r at all radii, and therefore 
a/af p = 1 and acceleration would be inefficient. 

In analogy with the previous discussion for equatorial field 
lines, let us now imagine uniformly expanding or contracting 
the field lines near the axis. That is, for each field line with 
a given 9f v <C 1 , let the polar angle far from the star become 
9 = k9f p , with the same value of k for all lines. For such a 
uniform expansion or contraction of the field, B p transforms 
to B p /k 2 . However, at the same time R becomes kR, and so 
a = B p R 2 is unaffected. In other words, there is no effect on 
acceleration. 

The key to obtaining acceleration along polar field lines 
is not uniform lateral expansion (divergence) or contraction 
(collimation) of field lines, but differential bunching of field 
lines. To see this rewrite equation (TT~9b as 



2 



1 



a 



1 



itB p R 1 



(26) 



where we have used the fact that af p « <&/-k near the pole. 
Clearly, for efficient acceleration, we must make B p substan- 
tially smaller than the mean enclosed field $ / nR 2 . That is, the 
field lines interior to the reference field line must be bunched 
in such a way that most of the flux has been pulled inward. 
As an example, consider a power-law distribution of the field 
strength, 

B p (R)<xR-t, (27) 



where the index £ measures the degree of bunching. This dis 
tribution gives 

/! 2 ' 



(28) 



which shows that the acceleration efficiency increases with 
increasing £. We reach equipartition between Poynting and 
matter energy flux (a ~ 1, 7 ~ /u/2) for £ = 1, and we obtain 
arbitrarily large efficiency (a <C 1) as £ — > 2. 



As we have described in §|3] polar field lines in simu- 
lation M90 are very efficient with -f/fi — > 1. According 
to the above, this would seem to suggest that £(/?) should 
be w 2 at the axis and should decrease with increasing 
R. The simulations, however, show that this does not hap- 
pen: as we show later, £ w 1, i.e., it stays roughly con- 
stant over a range of R. Instead, higher efficiency near 
the jet axis is achieved i n a different way: as field lines 
bunch around the jet axis dChiueh et al.l 1 1 99 1 : Eichlerj |T993t 
iBogovalovl 1 1 995b lBogovalov& Tsinganosl|1999h , they form 
a con centrated core (iHevvaerts & Normanlfl 989: Bogovalov 
1200 U iLyubarsky & Eichledl2001b iBeskin & Nokhrinal 120081) 
that takes up a finite amount of flux $0, leading roughly to the 
following poloidal field strength profile 



B p (R) = <f> 6(irR 2 )+B (R/R T 



(29) 



where we have approximated the core profile with the Dirac 
delta-function. This gives 



$o«ttV? 2 /(W/2), 



\i-7rzy?7$ , $ »^y?7(i-£/2). 



(30) 



That is, in the limit when the flux in the concentrated core is 
small compared to the flux in the surrounding power-law field 
distribution, the efficiency is the same as in ( 1281 ). However, in 
the opposite limit, i.e., sufficiently close to the core where the 
flux in the core dominates, the angular profile of the poloidal 
field distribution (and the value of £) become irrelevant for 
determining the acceleration efficiency. Thus, (1) differential 
bunching and (2) the resulting development of a concentrated 
core, are the key requirements for efficient acceleration. 

Note the following important corollary from the above dis- 
cussion. It does not matter whether the particular field line 
of interest collimates towards the axis or diverges from it. 
This has no effect on the acceleration. What we need is that 

(1) other field lines closer to the axis must converge more, 
or diverge less, compared to the reference field line, and/or 

(2) a concentrated core at the jet axis must contain a signifi- 
cant amount of magnetic flux. 

4.4. Comparison with Numerical Results 

Figure [5] shows results for models M90 and M10. The top 
panels show the behavior of a/ af p at different distances from 
the central star. We see that a / cif v decreases towards the pole, 
exactly where j/fj, is largest and a is smallest in Fig. [2] Also, 
Fig.|2]quantitatively confirms the validity of equation (fT9l by 
plotting prediction ( fT9l over the numerical solution. 

The lower panels in Fig.[5]illustrate the effects described in 
the previous two subsections. In the equatorial regions, we 
see that B p r 2 decreases from its footpoint value of unity. It is 
this decrease that allows whatever acceleration is observed in 
this region of the outflow. However, the decrease is modest, 
so the acceleration is not very large. 

For angles closer to the pole, B p r 2 actually increases rel- 
ative to its nominal initial value of unity. Nevertheless, this 
does not mean that there is deceleration because, as we argued 
above, acceleration near the axis is associated with differen- 
tial bunching, not with any overall expansion or contraction. 
For both M90 and M10, we see that the magnetic field sets up 
the required bunching so that the poloidal field is maximum at 
the axis and decreases with increasing distance from the pole. 
It is this outward decrease, coupled with the presence of mag- 
netic flux $0 in a concentrated core at the pole (see eq. [30b , 
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FIG. 5. — The panels on the left correspond to model M90 and those on the right to M10. The five curves correspond to different distances from the star: 
r = 10 2 , 10 3 , 10 4 , 10 s , 10 9 (see Fig.f2]for a description of line colors and types). The top panels show the quantity a/af„, which determines the acceleration 
along a field line. Note that ajai v is smallest near the poles, where Fig.f2]shows the largest acceleration. The middle panels show the profiles of normalized field 
strength B p r 2 at the same five distances. Notice the nonuni form distribution as a function of angle. This is necessary for acceleration, as explained in j|4.3l In 
the middle-left panel the curved dotted line corresponds to 13 It and accurately describes the angular profile of B p at large distances. The lower panels show the 
angular distribution of the flux function <£>/7r. 



that is associated with a decrease in a /dtp for polar field lines 
and the reason for strong acceleration. 

The effect is most clearly seen in model M90, as illustrated 
in Fig. Between r = 10 2 and r = 10 3 the angular profile of 
Bp becomes a steeper function of polar angle, leading to an 
associated decrease in a/flfp. Between r = 10 3 and r = 10 4 the 
trend reverses, and the poloidal field profile actually becomes 
a shallower function of polar angle; however, a / «f p continues 
to decrease. Beyond r ~ 10 4 the magnetic field profile does 
not evolve further. In this regime it can be well-fitted by a 
broken power-law, 

B p r 2 = 0.85 + 0.075 sin -11 6. (31) 

Therefore at these distances a = B p r 2 sin 2 9 does not evolve 
with r either. Despite this, a/ofy « 7ra/$ decreases with in- 
creasing r and the value of <£> increases. This is solely due to 
the increase in the amount of magnetic flux $o contained in 



the concentrated core. 

The picture that emerges is the following. The poloidal 
magnetic field establishes some equilibrium angular profile 
B p (r, 6)r 2 at intermediate latitudes that does not evolve with 
r. At progressively larger r, each individual magnetic field 
line becomes progressively more collimated towards the pole, 
but the same angular profile of magnetic field is maintained. 
This means that magnetic flux $ flows out from the low- 
latitude equatorial region (B p slightly decreases there), flows 
through intermediate latitudes without changing the profile of 
B p there, and ends up in the concentrated core at the polar 
axis, thereby uniformly shifting the angular poloidal flux dis- 
tribution $(#) up. The effect is clearly seen in the lower-left 
panel of Fig. [5] 

Note that it is not necessary to numerically resolve the con- 
centrated core in order to accurately describe the jet structure 
since it is only the total amount of magnetic flux contained 
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in the concentrated core that matters for the acceleration effi- 
ciency. In particular, our numerical method is well-suited for 
capturing such a core, even if unresolved, since our method 
conserves the magnetic flux to machine precision and can ac- 
curately capture the amount of magnetic flux that enters the 
core and remains there. We note that within a few grid cells 
from the polar axis, where the unresolved magnetic flux ac- 
cumulates, are least accurate but this does not affect the qual- 
ity of the solution at larger angles. To verify this, we have 
checked convergence of our models with angular resolution. 
For this, we ran a version of model M90 that uses a uniform 
angular grid and has a factor of 2 lower effective angular res- 
olution near the pole. Using this less-resolved model leads 
to a maximum relative difference in the flux function $ and 
Lorentz factor of less than 15%, even near the rotation axis. 
This difference is less than 2% at most radii (r < 10 2 and 
r > 10 5 ) and is smaller at larger 9. This confirms the accu- 
racy of the numerical solution. Future higher-resolution mod- 
els that resolve smaller angles and the concentrated core are 
required to determine how the solution connects the polar axis 
and for independent verification of results. 

5. ACCELERATION EFFICIENCY AND COMMUNICATION WITH 
THE AXIS 

5.1. Fast Magnetosonic Surface 

iBeskin et al.l (1998) have discussed the physical reason for 
inefficient acceleration in the equatorial regions. They show 
that it is related to the fast magnetosonic point. In the co- 
moving frame of a cold MHD plasm a, fast magnetosonic 
waves travel with a speed Vf given by (G ammie et alj|2003t 
McKinnev 2006b) 



lfVf = 



(32) 



where b is the comoving magnetic field strength. It is straight- 
forward to relate b to field components in the lab frame: 



,2 r2 f2 mb 2 mE\B v \m 

o — b — L £3 — - ?a — ^ pa. 



We thus find 



rV2 



(33) 



(34) 



Consider a streamline in the wind that moves outward 
with a local Lorentz factor 7. Let us first consider 
the limit of infi nitely high magnetization , p — > 0, i . e. the 
force-free limit dGol dreich & Juli an! 119691: [ Qkamotoi 1 1974 : 
Blandford 1976 : lLovelacdll976t iBlandford & Znaieklll977 
MacDonald & Thornd Il982t iFendt et al.l 119951: [K omissaro 
200 11 I2002allbh iMcKinnel l2006at iNarayan et al.l 120071 
TMNOl - In this limit 7/ — > 00, therefore throughout the 



solution we have 7 < 7/. This means that the fast magne- 
tosonic surface, defined by the condition 7 = 7/, where the 
wind becomes causally detached from fluid farther back along 
its streamline, is located at infinity. Such a force-free wind 
has a simple analytic solution in which the Lorentz f actor in- 
creases roughly linearly with distance (Michel 1973b, 



7: 



SIR. 



(35) 



In general, p^Q, yet we might expect that the behavior of 
the Lorentz factor in the sub-fast region 7 < 7/ is similar to 
the force-fr ee solution (f33T>. T his has been shown to indeed 
be the case (Beskin et al. 1998). However, the acceleration in 




lo gio r 



FIG. 6. — Dependence of various quantities as a function of distance for 
a field line with % p = it/ 2 in model M90. [Top panel] Solid lines show the 
radial dependence of 7, a, and fi from the simulation. [Bottom panel] The 
solid line shows the inverse of the fast Mach number 1 /Mf. The fast surface, 
indicated with the vertical dotted lin e, is located where Mf = 1 . This is also 
the location of the causality surface 14 U for thi s fie ld line. The various other 
dotted lines show the analytical approximation (49). 

the super-fast region 7 > 7f has been foun d to become log- 
arithmic, i.e. inefficient (Bes kin et al.| [T998). Once the wind 
has crossed the fast magnetosonic point 7 = 7/, it becomes 
causally detached from the fluid farther back along its stream- 
line. We might therefore expect efficient acceleration to cease 
beyond this fast magnetosonic point 9 for all field lines. 
If we define the fast Mach number Mf by 



M f = 



JfVf 



(36) 



then the fast point is the location at which My w 1 (this equal- 
ity would be exact if there was only motion along the poloidal 
field line; however, there is also a slow rotation in the toroidal 
direction which introduces a negligible correction that we ig- 
nore). For a relativistic flow (v 1), equation (fTTT) lets us 
recast the above expression in a useful form: 



2 3 
1 a (i— 7 



(37) 



Figure [6] shows results for a field line in M90 with 8f p = 6 = 
tt/2. Until the flow reaches the fast magnetosonic point, we 
see that a falls rapidly and 7 increases rapidly. However, both 
trends slow down substantially once the flow crosses the fast 
point. Beyond this point, Bes kin et al.l ((1998) have shown that 
a and 7 vary as the one-third power of log r. We confirm this 
dependence below. 

The relatively abrupt cessation of acceleration beyond the 
fast point for equatorial field lines is obvious in Fig.Q] where 
we see that 7 stops increasing once the flow crosses the fast 
magnetosonic surface (the middle of the three thick solid 

9 Note that fast magnetosonic waves move faster than Alfven waves, and so 
the causal horizon is determined by the fast waves rather than Alfven waves. 
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lines). However, it is also clear from Fig. [T]that something 
else operates on polar field lines. Model M10, in particular, 
shows substantial continued acceleration well after polar field 
lines have crossed the fast surface. We discuss next the rele- 
vant physics for these field lines. 

5.2. Communication with the Axis: Causality Surface 

We showed in §|4]that, for efficient acceleration along a field 
line, other neighboring field lines must shift laterally. At the 
equator, we need lines to move away from the midplane, while 
near the pole, we need field lines to experience differential 
bunching or develop a flux core. In order for any given field 
line to sustain efficient energy conversion, it must be able to 
communicate to other regions of the magnetized wind that un- 
dergo differential bunching or cause a concentrated flux core. 
This suggests that the fast magnetosonic point, which deter- 
mines where the fluid can no longer communicate back along 
its motion, is perhaps not so important. A more relevant is- 
sue is whether or not the fluid can communicate with regions 
near the axis that have field bunching or a flux core. We re- 
fer to the point at which a fluid element loses contact with 
the axis as the "causal point," and call the locus of causal 
points over all field lines as the "causality surface". By the 
above arguments, we expect that this surface, rather than the 
fast magnetosonic surface, plays the role of the boundary for 
efficient acceleration. 10 We expect efficient acceleration in- 
side this surface and inefficient, logarithmic acceleration out- 
side the surface (for a related d iscussion, see Zakamska et al. 
I2008t Komissa rov et al.l l2009t). In general fast waves prop- 
agate away from any given point toward the rotation axis 
through an intermediate region where the density, velocity, 
and magnetic field vary. Hence, one should trace the posi- 
tion of fast waves emitted from any given point outward over 
all angles and identify the "causal point" for each field line 
as where finally no such traces can reach the polar axis. For 
simplicity, we instead use only the local fast wave speed at a 
given point, and we identify the approximate "causal point" 
by where the locally emitted fast waves move away from any 
given point with a lab-frame local angle of 9 > with respect 
to the rotation axis, so that the waves do not reach the rota- 
tion axis over a finite propagation distance. We now calculate 
the approximate location of the causality surface using this 
approach. 

Consider a segment of the relativistic magnetized wind 
propagating with a velocity vector v and Lorentz factor 7 at 
an angle &j to the rotation axis. 11 In Appendix iDl we show 
that fast magnetosonic waves, which are emitted isotropically 
in the comoving frame of the fluid, in the lab frame will be 
collimated along v into a Mach cone with a half-opening an- 
gle 



municate with the axis, i.e., only if £ max > Oj, i.e., only if 



• * _ lf v f _ J_ 

71> Mf 



(38) 



By the argument given earlier, field line bunching and effi- 
cient acceleration are possible only when the fluid can com- 

10 We note that the causality surface is formally differ ent from the fast 
modifie d surface (which is discussed in de t ail in, e.g., | Guderle\ i 11962 ; 
IBlandford & PavnelfT982l ; IContop oulos 1995a; Tsinganos et al. 199jjj). Both 
of these surfaces are built on an idea of a full causal disconnect: the causality 
surface requires a causal disconnect across the flow, while the fast modified 
surface requires a casual disconnect along the flow. 

1 1 We use 9 for the polar coordinate of a point in the solution and 8 : for 
the angle between the local poloidal field and the axis. 



1 

w 



(39) 



For an equatorial wind, i.e., 9j = ir/2, equation (|39l shows 
that acceleration stops when Mf = 1, i.e., 7 = 7/ w <t 1//2 (c.f. 
eq. [34j assuming v/ — > 1). That is, the wind stops efficient 
acceleration as soon as it crosses the fast magnetosonic point. 
However, for smaller values of 9j, we obtain a different result. 

According to equation d38l . communication with the axis 
and efficient acceleration are possible until 



7« : 



lf v f 
sin 9 



rV2 



sint 



(40) 



The presence of the factor sin f9y in the denominator means 
that communication extends to larger values of 7, i.e., accel- 
eration efficiency becomes larger as 9j decreases. In other 
words, polar field lines can accelerate more easily. Using the 
definition of the fast magnetosonic Mach number d36i l. we ob- 
tain the following relation for the causality surface: 



sin 



{40} 1 03 <J C /"-7c 



' J.C 
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(41) 



where the subscript "c" indicates quantities evaluated at the 
causality surface. Inside the causality surface we expect the 
Lorentz factor to increase roughly linearly with distance, c.f. 
eq. d35l l. Based upon our earlier arguments, once outside the 
causality surface the acceleration will only be logarithmic. 
Using ( BTt . we can estimate the distance at which the causal- 
ity surface is located: 



1/3 



f2sin 5/3 e u 



(42) 



where we have assumed that (133b and 7 <C pL hold for r < r c . 

Figure [1] confirms that the causality surface (the outermost 
of the three thick solid lines) provides a better approximation 
to the boundary between the acceleration and coasting zones 
compared to the fast magnetosonic surface. Figure Q shows 
detailed results for a polar field line with sin f9f p = 0.1 in model 
M10. Notice that efficient acceleration continues well past the 
fast magnetosonic point ("F", dotted line on the left); accel- 
eration slows down only after the field line has crossed the 
causal point ("C", dotted line on the right). This is the reason 
why this particular field line is able to achieve a Lorentz fac- 
tor of 170 with a high efficiency of « 0.6, which is much 
larger than for equatorial field lines. 

5.3. Analytical Approximation 

We now develop an analytical approximation to calculate 
the Lorentz factor as a function of distance for any field line 
in a monopole magnetized wind. Generally, for a force-free 
jet , there ex ist two distinct acceleration regimes, as explained 
in TMN08. In Appendix [Al we generalize these results to 
MHD (finite-magnetization) jets. We summarize the results 
here. In the first acceleration regime, which is realized near 
the compact object, the Lorentz factor of the flow increases 
roughly linearly with distance: 



7^ 7o 2 + (^) 2 , 
where 70 is the initial Lorentz factor at r = 1 . 



(43) 
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FIG. 7. — Dependence of various quantities as a function of distance for a 
field line with sinf3f p = 0.1 in model M10. [Top panel] Solid lines show the 
radial dependence of 7, cr, and fi. [Bottom panel] The solid line shows the 
inverse of the fast Mach number 1/My and the dashed line shows the local 
opening angle of the field line Qi. The two vertical dotted lines indicate the 
positions of the fast magnetosonic point F (Mf = 1) and the causality point C 
(smdj = i/Mf) . T he various other dotted lines correspond to the analytical 
approximation 449t . 



Based on earlier arguments, beyond the causality surface 
the acceleration is only logarithmic. This is the second ac- 
celeration regime in which the Lorentz f actor is determined 
by the p oloidal shape of the field lines (B eskin et al.l 11998b 



TMN0 8|). Using the results of AppendixlBl as a guide (see also 
Beskin et al.ll 19981: iLvubarskv & EichleiEoOll) . we expect in 



this region 



72 oc In 1 / 3 r. 



(44) 



To make this formula quantitative, we demand that it gives the 
correct value of the Lorentz factor at the causality point r = r c 
(see eq.HTT): 

We do this by choosing the solution in the following form: 

l2 ~C nc ln l/3 (l+C 2 r/r c ) (45) 

- 1/3 ln 1/3 (l+C 2 r/r c ) 
~Cifl' . on „ ) ( 46 ) 



sin 2 / 3 ! 



J,c 



where C\ and C 2 are numerical factors of order unity that we 
later determine by fitting to the numerical solution. The radial 
and angular s calings in equ ation (l46l > a gree with the analytic 
expectation s (Beski n et al.lll998t ILvubarskv & Eichleril2001t 
lLvubarskvl2009t) . 

Note, however, that formulae (I44l -(|46]> become inconsis- 
tent at low magnetization since 7 cannot exceed [i, whereas 
the right-hand sides of these equations are unbound. Not- 



ing 



that the fast wave Mach number Mf is unbound and 



7/^i 1 / 3 k,M 2 J 3 for7 <C fi (eq.|37]i, we empirically modify d46ii 



FIG. 8. — Dependence of various quantities as a function of distance for 
a field line with sin6>f„ = 0.1 in model W10. Comparison to Fig. [7] shows 
that the wall inhibits the collimation of field lines. This brings the causality 
surface closer to the central star than in the model without a wall, and hence 
leads to a lower efficiency here compared to M10. [Top panel] Solid lines 
show the radial dependence of 7, cr, and fj,. [Bottom panel] The solid line 
shows the inverse of the fast Mach number 1 /Mf and the dashed line shows 
the local opening angle of the field line 8 j . The two vertical dotted lines indi- 
cate the positions of the fast magnetosonic point F (Mr = 1) and the causality 
point C (s'mBj = l/Mf). Th e various other dotted lines correspond to the 
analytical approximation )49t . 



by replacing 72//i 1//3 with M 2 J 3 - 



V 

2/3 ln I/3 (l +C 2 r/r c ) 
M7 =Ci- ' 



sin 



2/3 q 



(47) 



where C\ and C 2 are numerical factors of order unity (see be- 
low). Substituting for the fast wave Mach number Mf using 
equation (l37l i. we obtain a cubic equation for the Lorentz fac- 
tor on the field line beyond the causality surface, r> r c : 



1/3 ln 1 / 3 (l +C 2 r/r c ) 



sm 



2/3 q 



(48) 



where /1 is the value of the total specific energy flux on the 
field line in question (eq. 9j x is the value of the angle 6j 
that the field lines makes with the polar axis at the causality 
surface (see footnote ITTb. and C\ ~ C 2 ~ 1 (see below). In the 
limit 72 <C n this equation reduces to d46b . 

We now combine the two approximations (03]) and (08]) to 
write (see AppendixlAl 



1 



1 1 



T it 1 2 



(49) 



Clearly, the smaller of 7! and 72 determines the total Lorentz 
factor: in accordance with the above discussion, near the com- 
pact object 7 ~ 71 and at a large distance (outside the causal- 
ity surface) 7 « 72. We find that we obtain good agreement 
with our simulation results when we choose C\ = 2, C 2 = 0.4. 
In fact, for this single set of parameters formula (09]), with 71 
and 72 given by (03) and (08]), does quite well for all field lines 
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FIG. 9. — Dependence of various quantities as a function of distance for a 
field line with sin(9f p = 0.06 in model W5. [Top panel] Solid lines show the 
radial dependence of 7, cr, and fi. [Bottom panel] The solid line shows the 
inverse of the fast Mach number 1/My and the dashed line shows the local 
opening angle of the field line Qi. The two vertical dotted lines indicate the 
positions of the fast magnetosonic point F (Mf = 1) and the causality point C 
(smdj = i/Mf) . T he various other dotted lines correspond to the analytical 
approximation 449t . 

in all simulations, both in the limit of low and high magneti- 
zations. The various dotted lines in Figs. loT-flOl have all been 
calculated using this formula, and clearly provide an excellent 
representation of the numerical results. 

5.4. Other Models 

We have so far discussed in detail the representative models 
M90 and M10. However, we have carried out a number of 
other simulations. We mentioned models M20 and M45 in 
33.21 We have also carried out models with walls: W5, W10. 

Simulation W10 has the same setup as M10 but has an im- 
penetrable perfectly conducting wall at 9 = 10°. Figure [8] 
shows that the wall in simulation W10 keeps field lines near 
the wall from collapsing onto the pole and prevents the rest 
of the field lines from developing the lateral nonuniformity 
required for efficient acceleration. As a result, field lines 
very close to the wall accelerate more efficiently but the rest 
of the field lines have a suppressed efficiency: a field line 
with sin^fp = 0. 1 in model W10 has a lower Lorentz factor at 
r = 10 10 , 7 150, than the corresponding field line in model 
M10, which has 7 « 170 (compare Figs. |7]and[8]l. Simula- 
tion W5, which has the wall at 9 = 5°, is the most collimated 
model that we have simulated. Figure [9] shows a field line 
for that model that makes an angle sin6V p = 0.06 at the sur- 
face of the central star. This field line reaches equipartition 
by r ~ 10 I() with 7 w 200, the largest Lorentz factor we have 
achieved among all simulations in this paper. 

We have also performed a simulation called 1000M90 that 
has a uniform density profile on the stellar surface p(9f p ) = po 
with the same density at 9f p = as model M45. Thus, mod- 
els 1000M90 and M45 have similar profiles of density and 
fjb near the pole, but they differ near the equator. We expect 



the two simulations to show nearly identical behavior for po- 
lar field lines. This is indeed confirmed, as seen in Fig. [10] 
The point of this model is to verify that models with variable 
density, e.g., equation (fT2l . give reliable results near the axis, 
independent of how we modify the mass-loading of equatorial 
field lines. The numerical results confirm that this is indeed 
so. 

Figure [TO] shows transversal cuts through each of our mod- 
els at distances of 10 6 and 10 9 , and compares the simul ation 
results with the analytic approximation d49l described in 35.31 
For models with a constant density profile p(9f p ) on the star, 
the analytic approximation works extremely well at all dis- 
tances r and all polar angles. For models with a variable angu- 
lar profile of density on the surface of the star, the agreement 
is excellent along the field lines originating in the constant- 
density core while for other field lines, the best-fit values of 
factor C\ in equation d48| > apparently varies from one field 
line to the next. Therefore, adopting a single value C\ = 2 
provides only a rough description of acceleration along these 
field lines. 

6. DISCUSSION 

A number of previous authors have noted that efficient ac- 
celeration of a cold MHD wind requires field lines to diverge 
away from the equatorial plane. This field geometry was iden- 
tified as a "magnetic nozzle" due to the geometric similarity 
of jet nozzles intended to launch a supersonic flow. 

We generalized the concept of the "magnetic nozzle" by 
showing that the geometric bunching of field lines generally 
induces efficient conversion of magnetic to kinetic energy in 
neighboring regions. We clarified how this generalized "mag- 
netic nozzle" operates differently for equatorial and polar field 
lines. Near the midplane, lines merely have to diverge uni- 
formly away from the equator. The more they diverge, the 
larger the acceleration. For polar field lines, however, what 
is needed is neither simple divergence nor convergence, but 
differential bunching or a region that accumulates flux. A 
particular field line may either converge or diverge relative 
to its initial (purely radial) configuration. This has no effect 
on acceleration along this line. However, if neighboring field 
lines move such that the field strength decreases away from 
the rotation axis, e.g., as per the simple prescription given in 
equation d29l ), then acceleration will occur along the reference 
field line. The more the differential bunching (i.e., the larger 
the value of £ or $0), the larger the acceleration. 

The rearrangement of field lines described above requires 
different regions of the magnetosphere to communicate with 
one another, which is possible only if the flow speed is not 
too large. This introduces the second major difference be- 
tween equatorial and polar field lines. Equatorial field lines 
lose communication once their flow velocities cross the fast 
magnetosonic speed. This happens when the Lorentz factor 
7 ~ /i 1 / 3 , where p is the conserved energy flux per unit mass 
flux along the line. Equivalently, 7 ~ a 1 1 1 , where a is the lo- 
cal magnetization parameter. For highly relativistic flows, p 
is very large (e.g., p ~ 10 6 for the Crab Pulsar), so 7 is much 
less than p at the fast magnetosonic transition. Correspond- 
ingly, a ~ /i 2 / 3 is very large, which means that most of the 
energy flux is still carried as Poynting flux rather than as mass 
energy flux. 

For equatorial field lines beyond the fast point, a small 
amount of further acceleration is p ossible, but this on ly gives 
an additional logarithmic factor (Besk inet alj|1998h . After 
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FIG. 10. — Dependence of 7 (upper panels) and cr (lower panels) as a function of sin0f„ at r = 10 (left panels) and r = 10 (right panels) for a series of models. 
The names of the models are indicated next to the corresponding curves. Dotted lines show the analytic solution d49t and other lines show the simulation results. 
For models with a constant density profile p(8f„) on the star (M90, 1000M90), the analytic approximation works very well at all distances r and all polar angles. 
For models with a variable angular profile of density on the surface of the star (M45, M20, M10), the agreement is excellent alon g the field lines for which the 
initial density is constant; for other field lines, the agreement is less perfect, indicating the best-fit values of the factor C\ in equation d48t are different for different 
field lines. For models with the wall (W10, W5), the agreement is good near the axis and poor near the wall. 



allowing for this factor, the final asymptotic Lorentz factor 
on an equatorial field line at a large distance from the star is 
only 7 asyra - 2^ 1 / 3 ln 1/3 (r2r//i 1 / 3 ) (see eqs.|42]and|46). This 
is far smaller than the maximum Lorentz factor one would ob- 
tain if we had efficient acceleration along the field line, viz., 
7max = At- Thus, we confirm the previously known result that 
equatorial field lines in monopole geometry suffer from a se- 
rious a problem. We do not yet see any way of avoiding this 
conclusion. 

The situation is different for polar field lines. Even beyond 
the fast magnetosonic point, the fluid on these field lines can 
maintain communication with the axis (where field bunching 
allows efficient energy conversion). In fact, communication is 
lost only when 7 ~ cr 1 ' 2 / sinfy, where 9j is the angle between 
the poloidal component of the magnetic field and the rotation 



axis. For small values of 



0.01-0.1, as might be ap- 



propriate for relativistic jets, this gives a large increase in the 
asymptotic Lorentz factor reached by the flow. Including the 
additional gain from the logarithmic factor, we estimate (c.f. 
eqs.|42]and|48]l 

^, ^ Ctln^fflrsin^Z/xV 3 ) 

7asym«(M-7asym) ' T7T~~ , (50) 

- - / 



sin 



or, using 



CCD, 



Tasym 



l2 c\ l2 Wl 2 (nrsm 5 "ej, c /^) 
sin 9 1 r 



(51) 



where we note that [i — 7 asym oc a oc sin 9j C near the rotation 
axis. In the limit cr > 1 we can simplify d50l by approxi- 
mating /i-7as y m ~ These expressions are applicable be- 
yond the causality surface as approximately given by equa- 



tion d42l) . As an order of magnitude estimate, in these formu- 
lae one could substitute 9f p , 9, or 9j in place of 9 J)C . In any 



case, the sin 9j X factors in the denominator of equations ( pil l 
and ( 150b indicate that 7 aS ym is larger near the poles. 

6.1. Collimation and Acceleration 

Equations d40l >, (IBTT i predict that at roughly the same value 
of a ~ 1, an asymptotically more collimated simulation 
reaches a larger Lorentz factor: 



1 



'Tasym 1 



(52) 



where Ojr ~ 9f p > 9j is the angle at the causality surface. 
Figure QO] confirms this for the sequence of models M45- 
M20-M10, all of which have the same maximum value of fi, 
Mmax = 460. If all of these models could convert all of the elec- 
tromagnetic energy flux into kinetic energy flux, each of them 
would reach the maximum energetically allowed Lorentz fac- 
tor, 7 = /i max . While such a full conversion does not happen 
in any of these models, more collimated models reach higher 
Lorentz factors, in agreement with equation (l52l : each sub- 
sequent model in the sequence of models M45-M20-M10 is 
roughly twice as collimated as the previous one and at r = 10 9 
attains approximately twice as large a Lorentz factor, 50-95- 
170. 

6.2. Application to Relativistic Jets and GRBs 

The net conclusion of the previous discussion is that, 
whereas there is indeed a serious cr problem for pulsar winds, 
there is no similar problem for relativistic jets. We show that 
even if a flow is unconfined and has a monopolar-like shape, 
it can still efficiently accelerate in the polar region. For a jet 
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angle 6j ~ 2°, for instance, the scaling (IBlT l predicts an asymp- 
totic Lorentz factor 7 asym > 100 if the jet efficiently converts 
electromagnetic to kinetic energy and reaches a ~ 1 . In fact, if 
the jet is not efficient and carries more of its energy as Poynt- 
ing flux (e.g. for a larger value of fj,, see eqs. |5TT - l50l >, then 
a > 1, and we will have even larger values of 7 asym . These es- 
timates are confirmed by the numerical simulations described 
in this paper. 

The inferred total power of long GRBs is on the order 
of 10 51 erg dPiradl2005t fMeszarosll2006t iLiang et al1l2008l) : 
however, much less energetic events, with tot al energy release 
as low as 10 48 erg, have also been observed (So derberg et alJ 
2004, 2006). Let us compute the power output of jets in our 
numerical models and make sure that our jets are energetic 
enough to be consistent with these observations. The total 
power coming out from the compact object surface within the 
low-density core, 9 < 9 m . dx , is 

P jet = J 2Trr 2 d6sm6S(r,8)\ r= i m — J d9sm 3 6 



o 

n 2 4(2 + cos6» max ) . 4 







fl 2 

max a L riA 

<sin ~t~t " iax ' ( } 



where the last equality is for # max <C 1 . Converting the result 
to physical units, we obtain 



(54) 



where B r is the value of radial magnetic field component on 
the surface of the compact object and ro is the compact object 
radius. Evaluating the jet power for a magnetar with a charac- 
teristic period P = 1 ms and a surface magnetic field of 10 15 G, 
we get 



ni et , 



j2 x 10 4 
lms 



erg 



10 15 G 



ro 



10km 



u max 



(55) 



For a maximally-spinning black hole with dimensionless spin 
parameter a = 1, mass 3M Q , and surface magnetic field 
strength 10 16 G (lMcKinnevll2005l) . we have 



pj et 



:5xl0 48 fl 2 



/ B r \ 2 / M \ 2 /0 max 
V10 16 Gy \3M Q ) \ 10° 



erg 



(56) 

Given a characteristic duration of 10-100 seconds for a long 
GRB, our numerical jets provide, for the black hole case, 
10 50 -10 51 erg per event for model M10, 10 51 -10 52 erg for 
M20, and 10 52 -10 53 erg for M45. Therefore, the simulated 
jets from black holes in models M10 and M20 are ener- 
getic enough and move at sufficiently high Lorentz factors 
(7 ^ 100) t0 account for most long GRBs and can certainly 
account for low luminosity events. The energetics is also right 
for short GRBs: the simulated jets output 10 49 -10 SQ erg dur - 
ing a characteristic event duration of 1 second (Nakar 2007). 
For the magnetar case the energetics is lower: 10 47 — 10 48 — 
10 50 erg/s for a sequence of models M10-M20-M45. There- 
fore, the simulated jets from magnetars can account for less- 
luminous long GRB events and most short GRBs. 

While the energy fraction in the polar jet is small as com- 
pared to the total energy extracted by magnetic fields from the 



black hole, pJ et /p tot ~ 0^ (eq.[54] which assumes a uniform 
magnetic field distribution at the BH horizon), the absolute 
value of jet power PJ et is sufficiently large to account for long 
and short GRBs. We point out that the magnetic flux in ac- 
creting black hole systems is n on-uniformly con centrated in 
the polar region of the BH (e.g., McKinney 2005, due to am- 
bient pressure of the accretion flow), and therefore the total 
energy losses of the spinning black hole are actually domi- 
nated by the losses from the polar region rather than from the 
midplane region, meaning a larger fraction of power in the jet 
than given by eq. d54b . 

A very interesting question is whether the models suggest 
any characteristic value for the quantity jjOj: is this quan- 
tity generally smaller or larger than unity? For a jet with 
jjOj ^> 1 (say, ~ 10) only part of the jet within the beaming 
angle 6>/, w 1 /jj <c Oj is visible to a remote observer. As the 
interaction with the ambient medium decelerates the jet and 
the beaming angle becomes comparable to the jet opening an- 
gle, the edges of the jet come into sight and the light curve 
steepens achrom atically, displaying a "jet break" (Piran 2005; 
Meszaros 2006). In the other limit, jjOj <C 1, a jet would 
be incapable of producing achromatic breaks in GRB light 
curves. 

Achr omatic bre aks ha ve been fo u nd in pie-Swift 

times dFrail et al.1 l200ll: iBloom et al.1 120031: iZeh et alJ 
2006) and have yielded 8j ~ 0.03-1 rad. The situation 
is quite different in the post-Swift times for which a large 
amount of data available, and no jet breaks have been found t o 
fully satisfy closure relations in all bands (Liang et al. 2008). 
However, if one or more of the closure relations are relaxed, 
some of the breaks may be interpreted as "achromatic" and 
be used to derive jet opening angles tha t span a similar range 
as pie-Swift GRBs (Lian g etaLl 12008). Overall , it appears 
that the Lorentz factor of most GRBs is t,- > 100 (lPiranll2005t 
iMeszar os 2006, up to ~ 400, iLithwick & Saril 120011) which, 
with the above estimates for Oj, gives jjOj > 3 and indicates 
that in principle achromatic jet breaks are possible. 

Our numerical models have jj9j ~ 5-15, where the index j 
denotes quantities evaluated at the jet boundary which we de- 
fine as the boundary of matter-dominated region a < 1 . In this 
respect it is particularly fruitful to co mpare simulations M10 
and W10. As we discussed in ^5.41 the wall in model W10 
prevents field lines from collapsing onto the pole as much as 
they do in model M10. Comparison of field lines with a ~ 1 
for these models (see Figs. [7] and [8]), reveals a difference in 
the Lorentz factor of at most 20% and a much larger differ- 
ence in the value of jjO j (caused by a large difference in Oj 
due to the effect of the wall): "/j9j w 5 for M10 and « 12 for 
W10. The most collimated model W5 produces an even larger 
value, jjOj « 15. 

We now analytically confirm that in general the quantity 
jjOj ~ 5-10. For the jet boundary, using dBTI ) and charac- 



teristic values £1 



1. 



, 0.04-0.4, 



10' 



j )C ■ - U .UT- U .T, 

1 (see Figs. [7] [8] and|9]l, we get 



10 b -10 y 



7 A~ga>/ 2 log! ( / 2 ( 



10 2 -10 3 



(4-10)ct 



1/2 



(57) 



which is in good agreement with the simulation results. The 
large value of jjOj in this analysis arises solely due to the 
logarithmic factor that appears because a significant fraction 
of the acceleration occurs after crossing the causality surface, 
in the inefficient acceleration region: this is the case for all 
unconfined flows studied in our paper. This should be con- 
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trasted with confined outflows, collimated by walls with pre- 
scribed shapes, for which most of the acceleration tends to 
complete b efore crossing the causa lity surface and which have 
jjOj ~ 1 (Komissarov et al. 2009). We note that if we do not 
re quire that jets are matte r- domina ted, e.g. we allow a > 1 as 
in ILyutikov & Blandfordl 's (120031) model of GRBs, then the 
value of jj9j will be even higher (see eq.l57b. 

According to equation $5% , we can attribute the fact that 
some post-Swift GRBs show quasi-achromatic jet breaks, 
while many do not, by associating the former with jets that 
have jj9j > 1 and the latter with those that have jjffj < 1. 
Such a scatter in jjOj might be naturally produced by differ- 
ences in GRB environment (affecting Oj/0j, c > se e eq. |57Ti or 
the properties of the central engine (affecting a). Indeed, ac- 
cording to (1571 1. a low value of a (^C 1) or Qj (<§; 0j. c ) in our 
jets would mean jjOj < 1 and so the absence of a jet break. 

One could use equation ( f50b , which is based upon our ana- 
lytical model of the simulations, to obtain the Lorentz factor 
of any GRB jet. If the black hole or neutron star is nearly max- 
imally spinning with £1 ~ 0.25 and has a polar region with the 
reasonable value of /i ~ 1000 (TMN08), and if we consider an 
opening angle 6j ~ 4°, then by r ~ 10 s ~ 10 14 cm one obtains 
jj ~ 250 and jjOj ~ 17. For a range of opening angles with 
sufficient luminosity, one obtains a range of Lorentz f actors 
consi s tent with both s hort and long duration GRB jets dPiranl 
120051: lMeszarosll2006l) . Further, the product jjOj > 1, indi- 
cating an afterglow can exhibit the so-called "achromatic jet 
breaks," w here observations imply jj9j > 3 for long-duration 
GRB jets (lPiranll2005t rMeszarosl l2006). Our simulations and 
analytical models have jjOj ~ 5-15, which is proof of princi- 
ple that magnetically-driven jets can produce jet breaks. 

7. CONCLUSIONS 

We have studied relativistic magnetized winds from rapidly 
rotating compact objects endowed with a split-monopole 
magnetic field geometry. We used the relativistic MHD code, 
HARM, to simulate these outflows. We have constructed an- 
alytical approximations to our simulations that describe fairly 
accurately the Lorentz factor and the efficiency of magnetic 
energy to kinetic energy conversion in the outflow. 

Our main result is that, contrary to conventional ex- 
pectations, the winds from compact objects endowed with 
monopole magnetic fields have efficient conversion of mag- 
netic energy to kinetic energy near the rotation axes. We 
identify this polar wind as a jet since it contains a sufficiently 
high luminosity within the required opening angles of sev- 
eral degrees, and it accelerates to ultrarelativistic Lorentz fac- 
tors through an efficient conversion of magnetic energy to ki- 
netic energy (7 ~ n and a < 1 at large radii). We note that 
Lvubarsky & Eichler (2001) have identified a similar polar jet 
in unconfined magnetospheres based on its relative degree of 
collimation to the rest of the flow and the relativistic Lorentz 
factor. However, they did not concentrate on the accelera- 
tion efficiency and the transition to the matter-dominated flow. 
One can use equation d50b to show that, for example, order 
unity solar mass black holes or neutron stars with \i w a ~ 300 
near the compact object will readily produce jj ~ 150 at 10 14 
cm with 9j ~ 4° such that jjffj ~ 10. 

We are able to analytically explain how the jet efficiently 
converts magnetic energy to kinetic energy by identifying a 
"causality surface," beyond which the jet can no longer com- 
municate with the rotation axis that contains the flux core. 
When one region of the jet can no longer communicate to 
the flux core, that region ceases to accelerate efficiently. The 



communication between the jet body and the rotation axis al- 
lows magnetic flux surfaces and the Poynting flux associated 
with them to become less concentrated in the main body of 
the jet (at the expense of the bunch-up near the axis), and it 
is this process t hat a llows efficient conversion of magnetic to 
kinetic energy ( j|4.31 l. A similar mechanism, call e d "ma gnetic 
nozzle," was first described in Begelm an & Lil dl994l) . We 
clarify this mechanism by showing that the accumulation of 
flux near the rotation axis leads to a stronger decrease in a for 
polar field lines than for equatorial field lines. This effect is a 
ne w feature of idea l MHD winds that has not been discussed 
by Begelman & Li ( 1994) or any other a uthors except the very 
recent work bv lKomissarov et al.l (I2009I) . 

Our results demonstrate that ultrarelativistic jet production 
is a surprisingly robust process and probably requires less 
fine-tuning than previously thought. It is possible for spinning 
black holes and neutron stars to produce ultrarelativistic jets 
even without the presence of an ambient confining medium to 
collimate the jet. Further, we show that even unconfined (or 
weakly confined) winds from compact objects can produce 
sufficiently energetic jets {Lj ~ 10 49 erg/s) to explain many 
long GRBs and most short GRBs. 

We have confirmed the standard result that monopole mag- 
netospheres are inefficient accelerators in the equatorial re- 
gion. We have thus been unable to solve the er-problem for 
the Crab PWN under the assumption of an ideal MHD ax- 
isymmetric flow from a star endowed with a split-monopole 
magnetosphere. We note that at some distance from th e neu- 
tron s t ar an ideal MHD approximati on may break down (|Usovl 
Il994t ILyutikov & Blackmanl 120011) . No highly relativistic 
jet is observed in the Crab and Vela PWNe, despite our re- 
sults that suggest there should be such a feature. There are 
observations of non-relativistic jets with v/c ~ 0.5 that ap- 
pear diffuse and borderline stable. One way to resolve this 
discrepancy is that in PWNe systems the axisymmetric ul- 
trarelativistic jets we find are unstable to non-axisymmetric 
perturbations and so can be a prodigious source of high- 
energy particles and radiation via dissipation that causes the 
jet to slow to non-relativistic velocities (e.g., iGiannios et alj 
2009). T his notion of a visible jet emerging directly from 
the pulsar (Lvubarsky & Eichler 2001) is an alternative model 
to the more recent view that the observed jet is caused by 
a post-shock polar backflow with a > 0.01 that is guided 
by hoop stresses and force d to converge toward (and rise 
up along) the rotation axis (iKomissarov & L vubarskvl [2004t 
iDel Zanna et alj 120041) . In either model, one must con- 
sider non-axisymmetric instabilities since, in the backflow 
model, the shock structure and backflow could be highly non- 
axisymmetric and potentially unstable to non-axisymmetric 
instabilities. 

W e remark that wh i le pre paring this paper for publica- 
tion, [Komissarovetal] d2009) posted a paper describing ideal 
MHD simulations of confined and unconfined winds. They 
do make a minor note that their simulations of unconfined 
monopole outflows show efficient conversion of magnetic en- 
ergy to kinetic energy near the rotation axis. We find simi- 
lar results to theirs for the unconfined monopole wind. We 
are further able to make analytical estimates that explain the 
nature of this efficient conversion via introducing a causality 
surface at which the jet loses causal connection with the polar 
axis. We are also able to obtain a closed-form approxima- 
tion for the Lorentz factor in this region based upon a precise 
notion of the "causality surface." We note that previous au- 
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thors who studied stars endowed with a spli t-monopole field 
geometry in the ideal M HD approximation ( Bog ovalovl200 It 
iBucciantini et al.l [20061) did not perform simulations to large 
enough radii in order to observe the outflow achieving such 
an efficient conversion of magnetic to kinetic energy near the 
rotation axis or reaching such large Lorentz factors. 

We conclude with prospects for future research. Our so- 
lutions are axisymmetric, assume the compact object is en- 
dowed with a monopolar field geometry, assume particles are 
injected with 7 ~ 1 near the compact object, assume a 3> 1 
near the equatorial plane, and assume the ideal cold MHD ap- 
proximation holds. In future work we plan to consider the sta- 
bility of our solutio ns to non-axisymmetric instabilities (c.f. 
iNaravan et al.ll20 09), the effect of external confinement lead- 
ing to highly-collimated solutions, the injection of ultrarela- 
tivistic particles as suggested to occur, e.g., in the Crab pulsar, 
the injectio n of non-re lativistic particles a s may be relevant for 
thin disks (Bog ovalov & T singanosl fl999h . and the effect of 



a ~ near the equatorial plane as required near the pulsar's 
equatorial current sheet or for any system with an accretion 
disk. 
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APPENDIX 

A. LORENTZ FACTOR IN AN IDEAL MHD FLOW 



In this section we derive an approximate expressio n for the Lorentz factor in an id eal MHD flow. Conservation of energy and 
angular momentum flux along a field line imply (e.g.. lContopoulos & Ka zanas 2003): 



7(l-fLRiy) = 70, 

where is the angular frequency of field line and 70 and 7 are the initial and the local Lorentz factors of the field line, 

7 = (i-^-<r 1/2 . 



For convenience, introduce auxiliary variables x = flR and y: 



E=xB p , 
B v = -yB p . 



The drift Lorentz factor (TBI i is 



Resolving ((T), ( I A 1 b — d A4b for 7, we obtain: 



B 1 



l+y 2 



B 2_ E 2 1+y 2 



[ 7( i(iV)+*yf 



[l+y 2 -x 2 ) [xy(x 2 + 1 2 -\yl 2 + l0 {\+y 2 -x 2 ) l l 2 \ 



(Al) 

(A2) 

(A3) 
(A4) 

(A5) 
(A6) 



This equation is exact and expresses the MHD fluid velocity via the local parameters (x, y) and the footpoint parameters (70). 
In the limit x 3> 1, we have the following analytic approximation for the MHD Lorentz factor (assuming that \B V \ E, or that 



y 



2 -x 2 < x 2 ), 



lo+x 2 
1 +y 2 — x 



(A7) 



that is very close to the drift Lorentz factor ( IA5I ). This shows that asymptotically the Lorentz factor of the MHD flow ( IA6b 
is closely approximated by the drift Lorentz factor dA5] ) (c.f. lBeskin et al.| [l998. 2004: IVlahakisll2004 . We note that another 
convenient, albeit slightly less accurate, form of ( IA7b is 



% 2 : 



; 7o-l+7dr 



(A8) 



As we show by comparison to numerical results, formula ( IA71 > works extremely well not only asymptotically but at all distances 
from the star. We now show analytically why this is the case. In the monopole flow \BJ\ « E at all distances from the star (see 
34. Il l, therefore neglecting y 2 -x 2 as compared to x 2 , we obtain: 



j_ 

la 



7o 2 (l+x 2 )+x 4 



(j 2 +x 2 y/ 2 [7o+x 2 (7 2 +x 2 -l)V2] 



(A9) 



Accordi ng to this equation, the relative deviation of j a w.r.t. 7 is always smaller than 6% for any x = HR and 70. Therefore, 
formula (IA71 > can be used as an accurate approximation for the MHD fluid Lorentz factor at all distances from the star and initial 
Lorentz factor 70. 
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We now recast the approximation to the Lorentz factor 7,, iA7\ in a two-component form (c.f. lTMN08h : 



where 

increases roughly linearly with distance and 



-L«-L+-L (mo) 

7 7f 72 



1 \ = 1 l+ x 2 = 1 l+{mf (aid 



is related to the poloidal radius of curvature of the field line R c (TMN08) in the limit 72 -C /u: 

72~C^cos6» (A13) 

where C ~ 1 is a numerical factor. Since 7 2 /2 is the harmo nic m ean of 7^ and 7^, the smaller of 71 and 72 determines the total 
Lorentz facto r. Cl ose to the com pact object the first term in dAlOb dominates, and 7 increases roughly linearly with distance from 
the axis (see E2]and|™N0ll), 

7 « 7i = [7o + <mf] 1/2 « A/?. (A14) 

As we discuss in 35.21 beyond the causality surface the second term in (lAlOb becomes dominant and there the Lorentz factor 
grows logarithmically, 

7 » 72 oc ln 1/3 r. (A15) 
In the next section we present a compact derivation of this result. 

B. LOGARITHMIC ACCELERATION REGIME 

In this section we consider acceleration along near-midplane field lines, which have 9' = tt/2 — 9 <C 1, and show that it is 
logarithmic sufficiently far from the central star where the flow has switched to the second acceleration regime ( IA13b . We limit 
ourselves to a highly magnetized region of the flow, 7 <C /i. Assuming a near-monopolar structure in this region, we express the 
rate of change of angle 6' along a field line in terms of its local poloidal curvature radius R c as 



x 2 



d6' _ 1 JC9' 
dr R c rj 2 

where in the latter equality we used equation dA13b and the fact that r » R for 6' <C 1 . Differentiating equation d25l >. we express 
d r )Jd9 l along the field line: 

d9'~^~(9 J j 

Combining equations dBU and (1B2I ) and assuming a near-monopolar structure, 9' « we obtain a differential equation for 7: 



(Bl) 



M77^2- (B2) 



7 c?7 

- — -«C- , (B3) 
A* r 

or 

7«(C M ) 1/3 ln 1/3 (r/r»), (B4) 
where is an integration constant. In j j5.3| we argue that w r c , the radius at which the field line intersects the causality surface. 

C. AUXILIARY APPROXIMATE ANALYTIC SOLUTION 

HARM integrates the ideal cold MHD equations of motion by first interpolating the primitive quantities (i.e. density, velocity, 
and magnetic field) from cell centers to cell faces. T he cell faces will gen erally now have two values corresponding to an 
interpolation from different cell centers. As described in Gam mie"et al.l (12003), these two values are used to compute a generally 
dissipative flux that is used to advance the set of conserved quantities forward in time. These conserved quantities are then 
inverted to produce the new primitive quantities. 

In cases where one roughly knows the solution for the primitive quantities as functions of position, one can instead interpolate 
the ratio of the primitive quantities to the estimated solution. This reduces numerical dissipation by reducing the difference in the 
two values of primitive quantities interpolated to the face. The closer the estimated solution is to the true solution, the lower the 
numerical dissipation. So it is useful to estimate the solution as best one can. The particular form of this estimated solution has 
no other effect on the numerical simulation. 
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Here we provide our estimated solution based upon an approximate extension of known force-free solutions (TMN08): 



B r = r 2 , 


(CI) 


Be = 0, 


(C2) 


B^p = —QRB rj 


(C3) 


1 11 

7 2- 7 2 (OR) 2 ' ( M /2)2' 


(C4) 


2 2 i 
«r = 7 "I, 


(C5) 




(C6) 


r 


(C7) 


p=v—, 


(C8) 



with vectors given in an orthonormal (physical) basis and where r\ is given by the boundary conditions according to equation ©. 
We interpolate the numerical solution (given in a contravariant coordinate basis) divided by the estimated solution (also in a 
contravariant coordinate basis), with the excep tion of Be and ug for which we do not divide by the analytic solution. In addition, 
we indep endently interpolate 7 usin g equation ( 1C4| > in order to rescale the interpolated 4-velocity, which leads to a more accurate 
solution dTchek hovskoy et al. 2007)) • 

D. MACH CONE FOR FAST MAGNETOSONIC WAVES 

In this section we derive the expression for the half-opening angle of a Mach cone for fast magnetosonic waves. Consider a 
segment of the relativistic magnetized wind propagating with a velocity vector v and Lorentz factor 7. Let a fast magnetosonic 
wave travel in the comoving fluid frame at an angle £' to the direction v. Decomposing the wave velocity into parallel and 
perpendicular components, we have 

t>ll=«/COs£', = lysin^'. (Dl) 

Transforming back to the lab frame, the components of the wave velocity become 

VfCos£' + v Vfsm£' 

v\\ = — - , v ± = - . (D2) 

11 (1 + iitycos£') 7(1 +iw/cos4') 

Therefore, in the lab frame, the wave velocity vector is oriented at an angle £ to the vector v, where 

tan£=^= J ) (D3) 

vh 7(17- COS £' + 1>) 

In the lab frame, different waves have different values of £. Assuming that the medium moves with a superfast speed (i.e., 
v > v/), all waves move in the downstream direction in the lab frame and their £ values are all less than a certain maximum £ max . 
To determine £ max , we maximize tan£ with respect to variations of The maximum is achieved when 



v 2 — vh 



v 

This gives the half-opening angle of the Mach cone, £ m; 



Vf \ f 
cos£'= — sin£' = -i . (D4) 



tan£ max = — V * , (D5) 

sin£ max =^ = -|-. (D6) 

7« Mf 
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